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I.  INTRODUCTION 


In  the  analysis  of  stresses  in  gun  tubes,  particular  difficulties  are 
encountered  near  surface  boundaries  such  as  at  the  muzzle,  or  near 
discontinuities  in  stress  loading  such  as  at  the  rotating  band  of  the 
projectile.  The  analytic  description  of  stress  behavior  in  such  regions  is  in 
the  class  of  problems  termed  end  stresses  in  semi-infinite  cylinders. 
Currently,  the  Ballistic  Research  Laboratory  (BRL)  is  supporting  work  in  this 
area  through  contract*  with  Dr.  David  A.  Spence  of  Oxford  University  and  the 
Imperial  College,  U.K. 

The  research  discussed  in  this  report  was  conducted  in  support  of  Dr. 
Spence's  work.  Specifically,  in  view  of  the  existence  of  BRL'S  Bessel 
function  subroutine^  and  experience  in  computing  characteristic  roots  for  the 
hollow  cylinder,  Dr.  Spence  requested  the  BRL  to  provide  him  with  a  tabulation 
of  the  first  1000  roots  of  the  characteristic  equation  for  a  range  of  values 
of  Poisson's  ratio. 

II.  THE  CHARACTERISTIC  FUNCTION  FOR  THE  SOLID  ROD 

Cauchy's  method  of  characteristics  was  applied  by  V.K.  Prokopov^  to 
obtain  the  characteristic  equation  for  a  hollow  cylinder  with  stress-free 
surface.  A  similar  process  for  the  geometry  of  the  semi -infinite  solid  rod 
with  a  stress-free  surface  leads  to  the  characteristic  equation 


* European  Research  Office,  U.S.  Army ,  Contract  No.  DAJA  37-80-C-0192 . 

2K.L.  Zimmerman,  A.S.  Elder,  A.K.  Depue,  " User's  Manual  for  the  BRL  Subroutine 
to  Calculate  Bessel  Functions  of  Integral  Order  and  Complex  Argument,"  Tech¬ 
nical  Report  ARBRL-TR-02068,  Ballistic  Research  Laboratory,  1978  (AD  A056369). 

P 

V.K.  Prokopov,  "Equilibrium  of  an  Elastic  Axisymmetrically  Loaded  Thick- 
Walled  Cylinder,"  Prikadnaya  mathematika  i  mehanika,  Vol  XIII,  1949,  pages 
135-144.  Institute  of  Mechanics  of  the  Academy  of  Sciences,  USSR,  FTIO 
Translation  No.  J-2589 ,  Aberdeen  Proving  Ground ,  Maryland,  Translation  dated 
22  August  1967. 
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A(z)  «=  0 


(1) 


where 


zA(z) 


IQ(z)  [z2  +  (2-2v) ]  Ij(z) 

I^z)  z2  [lQ(z)j 


(2) 


Iq  and  1^  are  modified  Bessel  functions  of  the  first  kind,  of  order  0  and  1, 
respectively,  and  v  is  Poisson's  ratio. 

The  roots  of  Eq .  (1)  correspond  to  the  poles  of  the  integrands  in  the 

Fourier  integrals  representing  the  stresses.  The  values  of  z  satisfying  Eq . 

(1)  are  therefore  needed  for  calculations  based  on  residue  theory.  It  should 

be  noted  that  if  ZQ  is  one  such  root,  then  ZQ,  -Zq,  and  -ZQ  are  also  roots, 

where  Zq  denotes  the  complex  conjugate  of  Zq.  Thus,  only  roots  in  the  first 
quadrant  need  be  considered. 

To  evaluate  z  A(z),  use  was  made  of  the  BRL  Bessel  function  subroutine1  in 
the  equation 


zA(z)  =  z2[lQ(z)]2  -  (z2  +  C2~2v) ) [ Ix  (z2)]  .  (3) 

The  initial  estimate  used  for  the  distance  between  the  zeros  along  the 
imaginary  axis  was  tr  units.  Since  Newton's  method  worked  very  well  in  the 
simular  problem  for  the  hollow  cylinder,^  it  was  decided  that  this  method 
would  be  used,  along  with  the  initial  root  estimates  for  the  cylinder. 

III.  THE  FAILURE  (W  NEWTON'S  METHOD 

Once  Newton's  method  begins  to  converge,  it  converges  quite  rapidly.  It 
is  also  unfortunately  true  that  once  the  method  begins  to  diverge,  it  diverges 
rapidly.  This  fact  is  seen  at  once  from  the  expression  for  the  error 
terms.  We  have 


E 


k+1 


-f"(0) 

2f’(xk) 


<V 


(4) 


3R.W.  Hamming,  Numerical  Methods  for  Scientists  and  Engineers,  McGraw-Hill, 
Ina.,  Hew  York,  1973. 

4A.S.  Elder,  K.P.  Zimmerman ,  "Stresses  in  a  Hollow  Elastic  Cylinder  Produced 
by  a  Step  Function  of  Pressure  or  Shear 3 "  Memorandum  Report  BRL  MR  2454, 
Ballistic  Research  Laboratory,  1975  (AD  A008966). 
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where  is  the  error  in  the  current  estimate  x^  of  the  zero  of  the  function 
f,  is  the  error  of  the  next  estimate,  and  0  is  some  number  between  and 

the  true  zero  X.  If  f "( 0) /( 2f '(x^) )  is  nearly  constant  in  the  vicinity  of  the 
root,  then  the  error  at  each  approximation  decreases  (or  increases)  as  the 
square  of  the  previous  error. 

It  is  therefore  critical  that  the  initial  estimates  of  the  zeros  be  quite 
good  when  the  function  has  oscillations,  as  is  the  case  here.  In  the 
previously  cited  case  of  the  hollow  cylinder,^  at  most  9  iterations  on 
Newton's  method  were  required  to  obtain  each  of  the  first  five  hundred  roots, 
with  about  20-digit  accuracy.  The  characteristic  equation  for  the  hollow 
cylinder  contains  products  and  cross-products  of  the  four  modified  Bessel 
functions  Iq,  Ij ,  Kq ,  Kj ,  and  evidently  the  presence  of  the  K  function  tempers 
the  oscillations  somewhat.  In  addition,  a  rather  good  asymptotic 
approximation  to  the  characteristic  function  was  found,  providing  excellent 
initial  estimates  of  the  zeros. 

Such  was  not  the  case  for  the  present  work  with  the  solid  rod.  The 
initial  estimates  proved  to  be  insufficiently  accurate  to  provide  convergence 
in  Newton's  method. 


IV.  THE  MODIFIED  METHOD  OF  BISECTION 

3 

It  was  decided  that  a  root  search  based  on  the  method  of  bisection  might 
be  more  fruitful.  In  simple  terms,  the  method  of  bisection  consists  of 
dividing  the  region  near  a  zero  into  a  mesh  of  squares.  As  a  square  is 
traversed  in  the  counterclockwise  direction,  one  evaluates  the  function, 
seeking  a  change  in  sign  of  the  real  and/or  imaginary  parts.  This  signifies 
crossing  a  zero  curve  of  either  the  real  or  the  imaginary  parts,  or  both.  At 
this  point,  the  search  interval  is  halved,  and  the  process  continues  in  a 
direction  90°  counterclockwise  from  the  direction  just  used,  until  another 
sign  change  is  encountered. 

Because  of  the  analyticity  of  the  function  the  zeros  of  which  we  sought, 
we  did  not  need  to  concern  ourselves  with  most  of  the  caveats  about  the  use  of 
bisection  techniques.  The  location  of  the  first  several  zeros  was  estimated 
graphically.  This  was  accomplished  by  computing  functional  values  in  the  xy- 
plane  and  noting  sign  changes  in  the  real  and  imaginary  parts.  The  level 
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curves  U  =  0,  V  =  0,  where  U(x,y)  =  Re  {A(x  +  2y)}  and  V(*,y)  =  Im  {a  (x  + 

2y)  }  are  shown  in  Figure  1. 

From  this  point,  computation  proceeded  by  estimating  the  roots  to  be 

separated  by  it  units  in  the  direction  of  the  imaginary  axis.  Although 

somewhat  more  cumbersome  and  time-consuming  than  Newton's  method,  the  method 

“20 

of  bisections  provided  values  of  z  for  which  |a(z)|<1.0  x  10  on  the  BRL 
CDC  7600  computer.  Graphs  of  the  real  and  imaginary  surfaces  in  the  first  and 
fifth  quadrants  are  shown  in  Figures  2  and  3,  respectively.  Detail  was 
provided  by  alternating  color  contour  In  these  surface  plots. 

V.  CONCLUSIONS 


This  analysis  has  shown  the  method  of  bisections  to  be  viable 
in  the  method  of  characteristics.  Sensitivity  analysis  indicated 
change  beyond  the  first  10  roots  for  Poisson's  ratio  in  the  range 
This  information  was  forwarded  along  with  the  tabulation  of  roots 
Spence . 


tool  for  use 
little 
(0.1,  0.5). 
to  Prof. 
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Figure  3.  The  Imaginary  Part  of  A(z)  in  the  1st  and  5th  Quadrants 
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